Intermediate Image Analysis (optional)
Intro and Objectives
Hi all,
Here’re some more functions from the R stable for image analysis. Since I’m not sure if we’ll have time enough for this, am labeling this ‘optional’ only.
In addition to EBImage package we’d used, I’ll bring in OpenImageR as well this time.
The idea is to cover the following:
- Image filtering
- Image segmentation
There’s quite a bit more to cover but in the interest of the breadth vs depth tradeoff, I’ll stop here.
## Setup chunk
suppressPackageStartupMessages ({
if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") }
if (!require(EBImage)) { BiocManager::install("EBImage", version = "3.8")}
library(EBImage)
# if (!require(OpenImageR)){install.packages("OpenImageR")}
# library(OpenImageR)
})
# set this to your home folder
setwd(".")
getwd()## [1] "C:/Users/ANISHA/Downloads/lec_4_final"Filtering Images
Oftentimes, images have local artifacts and noise in them. To clean up this noise, we “smoothe” each pixel according to its surroundings.
An intuitive approach is to define a window of a selected size around each pixel and take the mean/median/some function of pixel values within that neighborhood.
After applying this filtering procedure to all pixels, the new, smoothed image is obtained.
EBImage wrapper function gblur() performs Gaussian image smoothing with filter size automatically adjusted to sigma for convenience. Behold.
require(EBImage)
# read in sample data
parrots = readImage("./parrots.jpg")
EBImage::display(parrots)
Above is the original image. Scroll your mouse over the image above and see what happens. Time now to apply the low-pass filter. See below.
## apply Gaussian low pass filter
parrots_gblur = gblur(parrots, sigma = 5) # 0.19 secs runtime
display(parrots_gblur, all=TRUE )
Notice how parrots_gblur appears blurred. In signal processing, such an image smoothing op is called low-pass filtering, i.e., only relatively low intensity pixel values are allowed to pass.
High-pass filtering is the opposite operation which allows to detect edges and sharpen images. This can be done, for instance, using a Laplacian filter (made explicit by defining the weights matrix around pixel).
Note: I’ll use EBImage’s linear 2D convolution filter - filter2() - to execute the high-pass filter. Try ?filter2.
P.S. We won’t go into what are convolutions and how the filters work, at this time. In DC these topics are only for exposure and illustrative purposes.
# defining laplacian filter's matrix kernel
lapl = matrix(1, nrow = 3, ncol = 3)
lapl[2, 2] = -8
lapl # view the laplace weights matrix## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 -8 1
## [3,] 1 1 1Think of the above as a weights matrix applied to each pixel and its surrounding 8 neighbours. The central cell in position c(2,2) is for the focal pixel.
Hence, in the ‘plain’ segment of an image (where neighboring pixels are similar to the focal pixel), the laplacian filter has the effect of cancelling out pixel values and rendering that segment dark.
It has the opposite effect at segment boundaries where pixel intensity changes sharply, suppose our focal pixel is at a segment boundary. Then the laplacian filter will cancel out pixel values on either side of the boundary while accentuating the intensity values of the boundary’s pixels.
P.S. Again, we don’t have to get into details at this stage but an intuitive understanding of what’s going on would be useful to develop at this stage.
Enough talk. Let’s apply the Laplacian filter and see for ourselves the results.
# using the FFT based convolution filter func 'filter2()'
img_lapl = filter2(parrots, lapl) # try '?filter2'
display(img_lapl)
Median filtering
Another approach to perform noise reduction is to apply a median filter, which has the advantage of removing noise while preserving edges.
The local median filter works by scanning the image pixel by pixel, replacing each pixel by the median on of its neighbors inside a window of specified size. This filtering technique is provided by the function medianFilter().
In what follows below, we first introduce random noise into the image –> then de-noise it using the median filter –> recover the ‘original’ image –> compare it with the original Image.
Behold.
## Inserting noise into the image.
n = length(parrots)/10 # n is a tenth of all the pixles in the image
Sampled_pixels = sample(length(parrots), n) # Drawing a random sample of 10% of the pixels
## now perturb the sampled pixels
parrots_noisy = parrots
parrots_noisy[Sampled_pixels] = runif(n, min=0, max=1) # runif() draws randomly from a uniform distbn
## view the result
display(parrots_noisy)
What does the pattern of the noise tell you? How can we smooth this now and revert to our original-ish image?
Let’s head there!
system.time({
parrots_median = medianFilter(parrots_noisy, 1)
EBImage::display(parrots_median)
}) # ~58 secs. computn takes time coz no FFT for nonlinear median computn.
## user system elapsed
## 64.72 0.04 64.98P.S. The above took time because median filtering is an inherently nonlinear technique and the FFT algo does not provide much help here. Again, no need to get into the details at this time.
Image Segmentation
‘Segmentation’ in Marketing refers to partitioning the customer base in to relatively homogeneous segments which share similar needs and preferences. Clustering methods are popular for segmentation such as the k-means algorithm.
Analogously, Image segmentation performs partitioning of an image, and is typically used to identify objects in an image. We’ll quickly see some illustrations next for demonstration and exposure sake, and move on.
# Using EBImage::bwlabel()
img = readImage("./india-night-lights.jpg") # india-night-lights
EBImage::display(img)
Above is a famous image. What can we glean from it?
# apply bwlabel
img_label = bwlabel(img) # 0.01 secs
max(img_label) # num of distinct objects found## [1] 3932# see output - integers corres to segments/ object groups
if (max(img_label) > 10){ table(img_label)[1:10]} else {table(img_label)}## img_label
## 0 1 2 3 4 5 6 7 8 9
## 512592 36 3 3 6 24 6 3 63 15bwlabel assigns positive integers as labels to each ‘segment’ it identifies. Table above shows how many pixles fall under which segment.
img_label # what's the obj like## Image
## colorMode : Color
## storage.mode : integer
## dim : 500 483 3
## frames.total : 3
## frames.render: 1
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0Let’s see what the segments are.
# display results
display(normalize(img_label)) # in grayscale
Above segmented image was n grayscale which can sometimes obscure segments.
We can view the above also in color by making each segment correspond to a different color. However since we have only a limited number of colors on the palette, the colors start repeating after a while.
See below.
display(colorLabels(img_label)) # in color## Only the first frame of the image stack is displayed.
## To display all frames use 'all = TRUE'.
Superpixel clustering
Next we use package OpenImageR to cluster larger or coarser segments comprising groups of similar-ish pixels into a ‘superpixel’ and plot the results.
See below.
require(OpenImageR)## Loading required package: OpenImageR## Warning: package 'OpenImageR' was built under R version 3.5.3##
## Attaching package: 'OpenImageR'## The following objects are masked from 'package:EBImage':
##
## readImage, writeImage#--------------
# "slic" method
#--------------
system.time({
res_slic = superpixels(input_image = parrots,
method = "slic",
superpixel = 200,
compactness = 20,
return_slic_data = TRUE,
return_labels = TRUE,
write_slic = "",
verbose = TRUE)
})## WARNING: The input data has values between 0.000000 and 1.000000. The image-data will be multiplied by the value: 255!
## The input image has the following dimensions: 850 567 3
## The 'slic' method will be utilized!
## The output image has the following dimensions: 850 567 3## user system elapsed
## 0.39 0.08 0.47str(res_slic)## List of 2
## $ slic_data: num [1:850, 1:567, 1:3] 176 178 179 180 182 186 190 193 196 198 ...
## $ labels : num [1:850, 1:567] 0 0 0 0 0 0 0 0 0 0 ...# display what we have now
plot_slic = OpenImageR::NormalizeObject(res_slic$slic_data)
plot_slic = grDevices::as.raster(plot_slic)
graphics::plot(plot_slic)
What do you see above? Notice how the partitions seem to have a similar size for the most part? And a similar composition within, for the most part?
Below is another func provided by OpenImageR package.
#---------------
# "slico" method [ the "slico" method does not require the "compactness" parameter ]
#---------------
res_slico = superpixels(input_image = parrots,
method = "slico",
superpixel = 100,
return_slic_data = TRUE,
return_labels = TRUE,
write_slic = "",
verbose = TRUE)## WARNING: The input data has values between 0.000000 and 1.000000. The image-data will be multiplied by the value: 255!
## The input image has the following dimensions: 850 567 3
## The 'slico' method will be utilized!
## The output image has the following dimensions: 850 567 3str(res_slico)## List of 2
## $ slic_data: num [1:850, 1:567, 1:3] 176 178 179 180 182 186 190 193 196 198 ...
## $ labels : num [1:850, 1:567] 0 0 0 0 0 0 0 0 0 0 ...And see the results below.
## display result
plot_slico = OpenImageR::NormalizeObject(res_slico$slic_data)
plot_slico = grDevices::as.raster(plot_slico)
graphics::plot(plot_slico)
Neat, eh?
In theory, there’s nothing stopping us from applying other clustering methods, obtaining clusters of similar pixels and grouping them together for display.
But all that is for a later time. Will stop here for now.
Sudhir